function MIo = cumPluvio(MI,nd,dt,d1,d2)   
%
%cumPluvio - Rainfall cumulative time series 
%
%   MIo = cumPluvio(MI,ND,DT)
%
% Let MI = [tI mIi] be a matrix such that mI(k) is the pluviometric output 
% value at the time tI(k), expressed in mm. 
% The output variable MIo = [t cum] is the cumulative rainfall at ND days
% computed with a step of DT hours, where the admitted DT values are 1, 4, 
% 6, 12 and 24, i.e. 1h, 4 h, 6 h, 12 h e 24h respectively. 
% If DT is undefined, empty or is not an admitted value, the default choice
% DT = 1 is used.
% The time series MIo is defined in the time span from ti = floor(t(1))
% to tf = ceil(t(end)). The values are provided at the times ti+n*DT, where 
% n = 0,1,2,...,(tf-ti)*24/DT. The n-th value are related to the range 
% time span [ti+n*DT-ND,ti+n*DT], boundaries included (this fact implies
% that MIo cannot be used as input data for other cumPluvio sessions).
%
%   MIo = cumPluvio(MI,ND,DT,d1,d2)
%
% In this case, MIo is computed in the time span from ti = floor(d1) to 
% tf = ceil(d2), where d1 e d2 are MATLAB serial date numbers.
%
% See also extraPluvio.

% G. Teza, 2020

if nargin < 5, d2 = []; end
if nargin < 4, d1 = []; end
if nargin < 3, dt = []; end

t = MI(:,1);
ti = t(1);
tf = t(end);

if isempty(d1)
    d1 = ti; 
end
if isempty(d2)
    d2 = tf; 
end
d1 = floor(d1);
d2 = ceil(d2);
if isempty(dt) || ~ismember(dt,[1 4 6 12 24])
    dt = 1;
    fprintf('Step: dt = 1 h\n');
end

na = 24/dt;     % numero di step nelle 24 ore
da = dt/24;     % step, espresso in giorni

to = d1:d2;
N1 = numel(to);
N2 = (N1-1)*na+1;               % N2 = N1+(N1-1)*(24/dt-1);
MIo = zeros(N2,2);
m = 0;                          % indice del file di output 
eps3 = 10^3*eps;                % per evitare "confronti duri"  
for k = 1:N1
    tk = to(k);
    if k == N1, na = 1; end     % per gestire il caso dell'ultimo elemento
    for kk = 1:na
        m = m+1; 
        tm = tk+(kk-1)*da;
        MIo(m,1) = tm;
        tm1 = tm-nd;
        tm2 = tm;
        Im = (t >= tm1-eps3)&(t <= tm2+eps3);
        Cm = MI(Im,2);
        Cm = nansum(Cm);    % Caso di pi� di un elemento
        MIo(m,2) = Cm;
    end
end           